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In this paper, we discuss a general procedure by which nonlinear power spectral densities (PSDs) of 
the harmonic oscillator can be calculated in both the quantum and classical regimes. We begin with 
an introduction of the damped and undamped classical harmonic oscillator, followed by an overview 
of the quantum mechanical description of this system. A brief review of both the classical and 
quantum autocorrelation functions (ACFs) and PSDs follow. We then introduce a general method 
by which the fcth-order PSD for the harmonic oscillator can be calculated, where k is any positive 
integer. This formulation is verified by first reproducing the known results for the k = 1 case of the 
linear PSD. It is then extended to calculate the second-order PSD, useful in the field of quantum 
measurement, corresponding to the k — 2 case of the generalized method. In this process, damping 
is included into each of the quantum linear and quadratic PSDs, producing realistic models for the 
PSDs found in experiment. These quantum PSDs are shown to obey the correspondence principle by 
matching with what was calculated for their classical counterparts in the high temperature, high-Q 
limit. Finally, we demonstrate that our results can be reproduced using the fluctuation-dissipation 
theorem, providing an independent check of our resultant PSDs. 

Keywords: power spectral density; harmonic oscillator; nonlinear optomechanics; quantum nondemolition 
measurement 


I. INTRODUCTION 


The harmonic oscillator, in which a particle is confined 
to a potential well that varies quadratically with position, 
has proven to be a very useful model in a number of clas¬ 
sical and quantum systems. In the classical regime, the 
harmonic oscillator provides an excellent description of 
periodic systems such as a mass on a spring or a pendu¬ 
lum, as well as resonating electronic LC circuits. In the 
realm of quantum mechanics, an analogous model is suc¬ 
cessful in predicting the behavior of a number of bosonic 
systems, such as photons confined to an optical cavity 
or phonons in an elastic solid. In fact, the vacuum itself 
is thought to consist of an array of harmonic oscillators 
with a broad range of frequencies pQ. 

Often, a harmonic oscillator model is applied to a sys¬ 
tem in isolation, where we generally consider only linear 
effects. However, when we begin to consider coupling 
between harmonic oscillators, or with other systems al¬ 
together, nonlinearities begin to enter the model, leading 
to new physics. An example of this sort of interaction 
arises in cavity optomechanics, in which two harmonic os¬ 
cillators, one describing an optical cavity and the other 
describing a mechanical resonator, are coupled to one 
another [2]. In this case, the motion of the mechanical 
resonator shifts the resonance frequency of the optical 
cavity, while the optics provide a radiation pressure force 
acting back on the mechanics. For moderate coupling, 
a simple linear model suffices, such that monitoring the 
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electromagnetic field provides a readout of the linear mo¬ 
tion of the oscillating mechanical device. However, as the 
interaction strength between the two systems increases, 
nonlinear coupling begins to occur, requiring that higher- 
order terms be added to the Hamiltonian 13®] • This 
provides a method by which one can obtain direct access 
to higher-order powers of the mechanical resonator’s mo¬ 
tion. For instance, a number of experiments have demon¬ 
strated direct coupling to the square of the oscillator’s 
displacement [3111EHII]. These types of measurements 
have generated significant interest, as they have been pro¬ 
posed as a method to perform quantum nondemolition 
(QND) measurements [121 S3 °f a mesoscopic quantum 
system [Sica El USEE!, as well as other exotic two-phonon 
processes, such as mechanical cooling/squeezing [S] and 
optomechanically induced transparency 


In order to make such measurements effectively, a 
knowledge of the autocorrelation functions (ACFs) and 
power spectral densities (PSDs) corresponding to the 
nonlinear readout of the oscillator’s motion is required. 
Though the first-order PSD is a well-known result pU- 
23, here we calculate a general PSD of any order for the 
quantum and classical harmonic oscillator, with a special 
focus on the linear and quadratic cases. The structure 
of this document is as follows. In Section [H] we provide 
a basic overview of the classical and quantum harmonic 
oscillators in the damped and undamped situations. Sec¬ 
tion III then provides a description of how to calculate 
the ACF and PSD of a classical time-dependent signal. 
Complementary definitions for a time-dependent quan¬ 
tum operator follow. Using the results of the previous 
two sections, Section IV introduces a general procedure 
that can be used to calculate the classical and quantum 
PSDs of fcth-order for the harmonic oscillator. Section 
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IVl reviews the k = 1 case of the first-order PSD of the 
harmonic oscillator, which is immediately followed by an 
extension to the k = 2 case of the quadratic PSD in Sec¬ 
tion [VlJ Finally, we conclude the document by discussing 
how these PSDs can be used in the context of real exper¬ 
iments. 


II. BACKGROUND 

A. Classical Undamped Harmonic Oscillator 

The model of the classical, undamped harmonic oscil¬ 
lator describes a system whose dynamics are governed by 
the following differential equation 

x + u>%x = 0, ( 1 ) 

where x(t) is a time-dependent variable that in this case 
we choose to be the position of the oscillator and ui 0 = 
yjkjm is the resonant angular frequency of the system, 
with k and m being the oscillator’s spring constant and 
mass, respectively. The familiar oscillatory solution to 
this second-order differential equation is given by 

x(t) = Xq cos(uJot + </>), (2) 

where xq and </> are an arbitrary amplitude and phase of 
the motion set by the initial conditions. 

We can determine the total energy of this system as the 
sum of its kinetic and potential energies. The potential 
for this system is V = \kx 2 = ^mco^x 2 , while the kinetic 

energy is simply the conventional K = | mx 2 = ,J ' m . 
where p is the linear momentum of the one-dimensional 
system. Using our solution for x(t) from above, we find 
the total energy to be 

I 

E = H = K + V= -h -?7 VjJqX 2 

2 to 2 

_ majoXo |^ s j n 2 + 4>) + cos 2 (wot + 4>)\ (3) 

rmjQX § 

= 2 ’ 

which is a time-independent quantity. Note that in this 
case, we can equate the total energy to the Hamiltonian 
of the system, which we have denoted as H. 


B. Classical Damped Harmonic Oscillator 

While the undamped harmonic oscillator provides the 
simplest solution to an oscillatory problem, this model 
can be made more realistic by introducing damping into 
the system, allowing for the description of real-world dis¬ 
sipative systems, including LRC circuits and nanome¬ 
chanical resonators [ 15 ] . The most straightforward way 


to introduce damping into Eq. ([I]) is to add a term pro¬ 
portional to x(t), producing the new differential equation 

x + Ti + ojqX = 0, (4) 


where T is a characteristic rate that quantifies the damp¬ 
ing in the system. 

In the underdamped case (T < 2w 0 ), the solution to 
this equation is given by 

r t 

x(t) = Xoe~^~ cos {ujdt + 4>), (5) 


where u>d = — (T/2wo) 2 is the shifted resonance 

frequency due to damping. When damping is very small 
(T -C 2wo), we neglect this shift and take ~ u>q. In this 
limit, Eq. ([2]) provides a good approximation for the solu¬ 
tion of the damped harmonic oscillator given by Eq. (|5j) . 
From this point forward, we will assume we are in the 
small damping limit, as this is the case of interest for 
most nanomechanical systems. 

Another useful parameter which can be used to quan¬ 
tify the damping of the system described above is the 
quality factor Q , defined by the equation 


q = 2t -Ke’ 


( 6 ) 


where A E is the energy dissipated per oscillation cycle. 
For the damped harmonic oscillator given above, we cal¬ 
culate the total energy of the system as we did in the 
undamped case, resulting in 


E = 


muinX, 


bUo -rt 


2^,2 


[sin 2 (w 0 t + </>)+ cos 2 (ix 0 t + </>)] 


( 7 ) 


o^o -rt 


where in the above equation we have neglected a term 
in x(t) that is proportional to T as we are in the small 
damping limit. The result for the energy of the damped 
harmonic oscillator is identical to that for its undamped 
counterpart given in Eq. ([3]), except now the energy de¬ 
cays on a timescale set by F. The energy dissipated in 
one cycle is then given by the change of energy over one 
period of oscillation tq, that is 


A E = ™ 2 o x2 o „-r t 


e — 


mulxl c _ rt 


(1 


mulxl „_r( t+rn ) 
2 

- e~ rT0 ) . 


( 8 ) 


The quality factor for this system is then given by 


Q 


2n 

1 - e~ rT ° 


27T LO 0 

1 - (1 - rr 0 ) = T’ 


( 9 ) 


where we have again used the small damping limit and 
the fact that we can relate the period of oscillation to 
the system’s angular resonant frequency via r 0 = 2tt/uiq. 
From Eq. it becomes apparent that the small damping 
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limit is equivalent to the high-Q limit, as smaller damping 
leads to a reduction in energy dissipation. In fact, using 
our condition for the small damping limit above, we can 
quantify the high-Q limit as Q 3> 1/2. For the remainder 
of the document, we will refer to the small damping limit 
as the high-Q limit. 

By analyzing the undriven, damped harmonic oscilla¬ 
tor above, we were able to investigate the time domain 
solution of the resonator’s motion, as well as its energy 
dissipation. However, this description is still somewhat 
incomplete as generally the motion will be driven by some 
time-dependent external driving force f(t). In such a sit¬ 
uation, we arrive at the driven differential equation of 
motion 


x + Ta; + ujgX = —. (10) 

m 

Analytical solutions for x(t) in this case can only be de¬ 
termined for a small number of special cases of /(f), such 
as a harmonic driving force. However, it is often more 
fruitful to Fourier transform this equation to get its ex¬ 
pression in the frequency domain, resulting in 


x(uj) = xM/(w) 


/ M 

m (wq — w 2 — *wr) ’ 


( 11 ) 


where x(w) = E{x(t)} and f(oj) = E{f(t)} are the 
Fourier transforms of x(t) and /(f) as defined in |A l| and 
we have used the property in Eq. ( |A4| ) to calculate the 
Fourier transforms of the derivatives. We have also in¬ 
troduced the generalized mechanical susceptibility 


xM 


1 

TO (UJ g — W 2 — ZwT) ’ 


( 12 ) 


which allows us to relate the resultant position to the 
applied force in the frequency domain. 


C. Quantum Harmonic Oscillator 

To extend the above treatment of the harmonic oscilla¬ 
tor into the quantum domain, we must first determine its 
governing quantum mechanical Hamiltonian. This is ac¬ 
complished by simply replacing x and p in the first line of 
Eq. ([3]) with the canonically conjugate position and mo¬ 
mentum operators x and p. With these new operators, 
our kinetic and potential energies now become K = J‘ m 
and V = ^uiojqX 2 , resulting in 


the eigenstates, however, the energies are given by 


E n = hio 0 



(14) 


where H = h/2n is the reduced Planck’s constant. In 
the above equation, n is an integer and signifies the state 
of the oscillator. Quantum mechanically, this number n 
can be interpreted as the number of quanta in the sys¬ 
tem, for example photons in a cavity or phonons in a 
solid. Therefore, E 0 = denotes the ground state en¬ 
ergy where n = 0 and no quanta exist in the system. 
This energy will be shared evenly between the expecta¬ 
tion values of the kinetic and potential energy such that 
(K) = (V) = ^ in the ground state. 

We now introduce the raising (creation) and lowering 
(annihiliation) operators E and 6, also known as the lad¬ 
der operators. These two quantities are given by 


b = 

Si = 


TOW 0 

~W 

rriLOo 

~2h 


x — 


mooo 

i 

—i 
mixo 


(15) 


and obey the commutation relation [6, E] = 1. These 
operators are convenient as they produce the following 
simple relations when operating on the energy eigenstates 
of the system 


b | n) = \fn\n — 1), 

Si | n) = Vn + 1 \n + 1), 

as well as their Hermitian conjugates 

(?r| b' = (n — 1| \fn, 

(n\b = (n + 1| \Jn + 1. 

From these relations we also have 
Eb |n) = n \n ), 
fofol | n) = {n + 1) | n). 


(16) 


(17) 


(18) 


As seen above, when acting on the energy eigenstates 
the operator N = Eb returns the number of quanta n 
of that state and is known as the number operator. By 
inspecting Ecp (14) it should therefore be clear that the 
Hamiltonian can be expressed as 


H = Hu>o 




(19) 


^ /s ^ rp 1 1 

H = K + V = --1- -toWqX 2 . (13) 

Inputting this Hamiltonian into the Schrodinger equa¬ 
tion, it is possible to solve for the energy eigenstates of 
this system |n), along with their corresponding energy 
eigenvalues E n , where H \n) = E n \n). For this deriva¬ 
tion, we do not concern ourselves with the exact form of 


We can also write x in terms of the ladder operators as 
x = x zp [ (b + , (20) 

where we have introduced x zp { = 2mui 0 1 wl^l 1 lh e 
amplitude of the quantum mechanical zero point fluctu¬ 
ations of the oscillator. 
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Up to this point, we have been dealing with opera¬ 
tors in the Schrodinger picture, where it is the eigen¬ 
states, not the operators, that carry the time-dependence 
of the problem. However, since we are dealing with time- 
dependent signals, it is convenient to turn to the Heisen¬ 
berg picture of quantum mechanics, where the operators 
are now the quantities that vary in time. The dynam¬ 
ics of an operator O, which was time-independent in the 
Schrodinger picture, is now governed by the differential 
equation 


o = - h [H M . 


( 21 ) 


Upon inspection of this equation, we see that an operator 
which is time-independent in the Schrodinger picture and 
commutes with the Hamiltonian will remain constant in 
the Heisenberg picture. 

Using the Hamiltonian for the quantum harmonic os¬ 
cillator given in Eq. (19), along with the commutation 


relation for the ladder operators, we can obtain a differ¬ 
ential equation for b(t) given by 


b = —icjob. (22) 

This equation is easily integrated to obtain the expres¬ 
sion for the annihilation and creation operators in the 
Heisenberg picture as 


bit) = be~ iuiot , 


(23) 


where the latter equation is obtained by taking the ad¬ 
joint of the former. 

Finally, we determine a time-varying expression for 
x(t) in the Heisenberg picture by inputting the relations 
in Eq. (23) into Eq. (20) to obtain 


x(t) = x zp{ (be~ iuot + . 


(24) 


Another advantage of working in the Heisenberg pic¬ 
ture is that we can introduce damping into Eq. (22) 
through a formalism known as input-output theory QJ . 
To do this, we assume that our harmonic oscillator is cou¬ 
pled to a bath which has some effective temperature T. 
Through this coupling, the oscillator is able to reach ther¬ 
mal equilibrium with the bath, by either losing energy to 
it or gaining energy from it, corresponding to damping 
of the oscillator and an incoherent drive from the bath. 
In general, this drive will have contributions originating 
from both the thermal occupation of the bath, as well as 
its quantum mechanical fluctuations. 

Often, this bath is chosen to be an ensemble of har¬ 
monic oscillators with varying resonance frequencies, all 
of which are at the bath temperature. This method 
proves to be very effective, as it is solvable due to the 
simplicity of the harmonic oscillator, and provides an 


accurate model of physically realizable baths, such an 
electromagnetic held or phonons in a solid [Tj. Using this 
model, along with the first Markov approximation (mem¬ 
oryless coupling to the bath quantified by a constant) [T], 
we modify Eq. (22) to obtain a new equation of motion 


S 7 = -icoob-y - —b-y + Vfb n , 


(25) 


where we have introduced a subscript 7 to differentiate 
this ladder operator from the undamped one. In the 
above equation, T quantifies the coupling of our oscil¬ 
lator to the bath and corresponds directly to the me¬ 
chanical damping rate mentioned above for the classical 
case. This is exemplified by the fact that if we set T = 0 
in Eq. (25), coupling to the bath is severed and we re¬ 


claim the original, undamped differential equation given 
by Eq. ( f22| ). 

In the above equation, the two new terms have arisen 
from coupling our oscillator to the bath. The second term 
on the RHS describes a decay in the amplitude of 6 7 (t) 
due to energy radiation to the bath, while the third term 
represents the drive due to input noise from the bath, 
given by the operator b n {t). We assume that this noise 
will be delta-correlated in time ( i.e. Markovian), which 
corresponds closely to classical white noise, resulting in 

mu 


( K(0 n (t ')) = (n b + 1) S(t - t'), 

= n. b 6{t - t'), 

(bUtPnit')) = (b n (t)b n (t ')) = 0. 


(26) 


Assuming that the bath occupation will be constant over 
the small bandwidth ~ T of interest about the oscilla¬ 
tor’s resonance frequency we can take the bath occupa¬ 
tion number to be the single value n b = n b (u) 0 ). 

In this case, it is difficult to obtain a time-domain rep¬ 
resentation for & 7 (f) due to the noise input into the sys¬ 
tem. Instead, we Fourier transform Eq. (25) to obtain 


the spectral form of the annihilation operator 


i r,.\ VTb n (uj) 

• ( \ -p /o 5 

1 (cjq — + 1/2 


(27) 


where we have introduced the Fourier transformed oper¬ 
ators & 7 (w) = .F{ 6 7 (f)} and b n (oj) = T{b n (i)}. We can 
also determine the spectral form of the creation opera¬ 
tor, b\(tjj) = .F{W)(f)}, by taking the adjoint of the above 
equation and using the relation [& 7 (w)]i’ = &j)(—w), which 
results in 


b\{ui) = 


VfbUoj) 


—i (wq + to) + T /2 


(28) 


Using these results for b 7 (u>) and b^ (w), we find a damped 
representation of the position operator to be 


x 7 (uj) = x zp[ ^b 7 (uj) + b\(u) 


(29) 
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Finally, with our definition of the inverse Fourier trans¬ 
form, along with Eq. (26 1 , we obtain the correlators for 
the Fourier transforms of the bath operators in frequency 
space as 


(b n (u)bl l (u')) = 27r (n b + 1 ) 8{lo + a/), 

(b'l(u)b n (u 1 )) = 27 m b S(uj + a/), (30) 

<&nM&n( w ')} = {bn(u)b n (u')) = 0. 

Note that a difference of a factor of 27 t arises between 
these correlators and others found in the literature HI 1231 
due to our definition of the Fourier transform. These op¬ 
erators with damping included will be useful later when 
determining the PSD for the damped quantum harmonic 
oscillator. 


III. AUTOCORRELATION FUNCTIONS AND 
POWER SPECTRAL DENSITIES 

In this section, we provide definitions that allow us to 
calculate the ACFs and PSDs for classical signals and 
quantum operators. Note that in this document, we will 
introduce a bar over the classical ACFs and PSDs to 
differentiate them from their quantum analogs. 


A. Classical 

We begin with a classical description of the ACF for a 
real, time-dependent signal a(t). The ACF tells us how 
the value of a{t) at a time t' is correlated to itself at a 
later time t + t' and is given by [25: 

1 r°° 

R aa (t) = lim — / a(t')a(t'+1) dt'. (31) 

T 0 ->oo T 0 J_ OQ 

Furthermore, by taking t = 0, that is inspecting how a(t) 
is related to itself at the same time, we obtain the time 
average of a 2 (t) defined as 

1 r°° 

(a 2 ) = Raa( o) = lim — a 2 {t')dt', (32) 

Tq— foo 1 0 J_ oa 


ACF from Eq. (31) to write this PSD in terms of the 


Fourier transform of a(t ) as 


S aa (w) = lim -^|a(w)| 2 , 

T 0 ->-oo 1 o 


(34) 


where a(w) = J-{a(t)}. By performing the inverse 
Fourier transform we can also recover the ACF from the 
PSD as 

f OO 

S aa (u)e- iut du. (35) 


1 


R aa {t) - ^ 


Also, through Eq. (32) it is apparent that the PSD is 


related to the time average of the squared signal by 

1 r™ - 

(a 2 ) = — J S aa (u)duj. (36) 

Generally, the energy of the signal is proportional to the 
signal itself squared, so by integrating the PSD over all 
frequencies, we are able to determine the average energy 
of the signal in question. This property will be useful 
later when normalizing our PSDs. 

We conclude our discussion on the classical PSD and 
ACF by noting that the definitions we have chosen are 
for the two-sided PSD, which is defined for both positive 
and negative frequencies. We elect to use the classical 
two-sided PSD for this document, as it is easier to cor¬ 
respond with the quantum PSD, in which an asymmetry 
between positive and negative frequency arises. How¬ 
ever, we mention briefly that when performing classical 
experiments, it is sometimes more convenient to work 
with the one-sided displacement PSD, which is defined 
over only positive frequencies [25] and is often quoted 
in the literature HEH 2 I]- Using the fact that a classical 
two-sided PSD is an even function, we can see that in 


performing the integrals in Eqs. (35) and (36), the limits 


can be changed from 0 to oo, provided we multiply by a 
factor of 2. Therefore, we can determine the one-sided 
PSD by multiplying the two-sided PSD by a factor of 
two and restricting its definition to be over only positive 
frequencies. This simple conversion from a two-sided to 
a one-sided PSD applies to all classical PSDs derived in 
this document. 


B. Quantum 


where we have used the shorthand (a 2 ) = (a 2 (t)) and 
will continue to use this notation throughout the docu¬ 
ment. 

The PSD, which specifies the signal’s intensity at a 
given frequency, and the ACF for a signal are related 
to each other by a Fourier transform. Therefore we can 
obtain the PSD for a(t) from its ACF by [25] 


In the realm of quantum mechanics, physical observ¬ 
ables correspond to Hermitian operators that act on 
wavefunctions. Therefore, our ACF and PSD will be in 
terms of the averages of these operators. 

The quantum PSD is a spectral function that tells us 
the intensity of a time-dependent quantum mechanical 
operator a{t) at a given frequency w and is defined as 


S aa (u)= R aa (t)e iut dt. (33) 

J —OO 

Furthermore, we can use the properties of the Fourier 
transform given in along with the definition of the 


S aa (w) = I R aa (t)e iut dt 

) 

(a(t)d(O)) e lult dt, 


( 37 ) 
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where R aa (t) = (a(t)a( 0)) is the ACF for a(t). At a finite 
temperature T, we can determine the ACF for a(t) from 


(a(t)a(O)) 


Tr {e~ pfl e iflt / n ae- ikt l n a} 
Tr {e~PH} 


(38) 


where H is the Hamiltonian of the system, /? = 1/fcgT 
with ks being the Boltzmann constant and Tr{} denotes 
the trace of an operator. In this paper, we choose to 
work in the energy eigenstate basis so that the trace of 
an operator O is given by 


Tr{6} = ]T (n\ 6 |n), (39) 

n 


where |n) is the nth energy eigenstate of our quantum 
system. We can therefore see that the denominator of 
Eq. |38|), given by 


-0E n 


z = ^2(n\e fiH |n) = ^ e (n|n) = ^ e 

n n n 

( 4 °) 

is the canonical partition function [I]. 

We also point out that we can inverse Fourier trans¬ 
form S aa (ui) to obtain R aa (t) as 

1 r°° 

Raa(t) = (o(t)o(O)) = — / S aa (uj)e duj. (41) 
Setting t = 0we then have 


(a 2 ) = ^ S aa {v) duj, 


(42) 


in direct correspondence with Eq. (36) for a classical sig¬ 
nal. 

It is also possible to express S aa (u>) in terms of the 
Fourier transform of a(t). By inputting the definitions 
for the Fourier transform of a(t) into Eq. (37) we find 


1 f°° 

Saa(v) = — J (a(w)a(u/)) du', 


(43) 


where a(w) = F{d(t)}. This relation is very useful, as it 
provides an alternate method by which we can calculate 
PSDs using the frequency domain. Equipped with these 
definitions, as well as the relations given in Section [TTJ we 
are now ready to determine the PSDs for the harmonic 
oscillator. 


IV. GENERAL FORMULATION FOR THE 
POWER SPECTRAL DENSITY OF x k 


A. Classical 


Beginning with the fcth-order classical PSD, we use the 
definition of the PSD as the Fourier transform of the ACF 
given in Eq. (33), along with Eq. (34), to obtain 


S x k x k(w) =J r {R x k x k(t)} = lim -* r |x (fc) (w)| 2 . (44) 

TO —>00 1 0 

Here we have used the notation x^^lo) = F\x k (t )\ t o 
denote the Fourier transform of x k (t). Using Eq. (A8), 
we can express this quantity as 

x( k \u) = F{x(t) ■ x(t) • ... • x(t)} 

1 r t t t ^ (45) 

= (■ 2n) k ~ l X ^ * X ^ 


where the ellipsis (...) indicates that the corresponding 
operation is performed on k terms (for a total of k — 
1 operations). We point out that with this notation, 
k = 1 corresponds to a single term with no operations 
performed. 

For a general driving force /(w), this expression is very 
difficult to solve. However, if we restrict ourselves to a 
frequency-independent drive ( i.e. f(u) = F), as is the 
case in thermally driven classical oscillators, the problem 
simplifies significantly, as we obtain the relation 

x(co) * x(u) * ... * x(uj) = F k x{u) * x( w ) * * X( w )i (46) 

where we have input the relation in Eq. for x(oj). 
We can then write 


S x k x k(uj ) = S% Fk |x(w) * xM * ... * xM| 2 , (47) 


where we have defined a white noise thermal force PSD 


S% Fk = lim — 
* To^oo T 0 


]?k 


(2tt) 


k-1 


(48) 


Note that while the driving force is constant in frequency 
space, it still grows as we increase T 0 , balancing out the 
division by infinity such that Fk . remains constant. 
The value of this quantity can be determined by ensuring 
that Eq. (36) is satisfied. In the high-Q limit, we can 
approximate the expectation value of (x 2kX ) for a damped 
harmonic oscillator as that for the undamped oscillator in 
equilibrium with a bath at temperature T (see[C]), which 
results in 


1 f°° 

(x 2k ) = ^ J S x k x k(u}) du> = x^ 


Ml 

2 k k\ ’ 


(49) 


We now introduce a method by which the PSD can be 
calculated in both the classical and quantum regimes for 
any power of the position of a harmonic oscillator x k (t ), 
where k is any positive integer. From this point forth, 
we label the PSD for x k {t) as the fcth-order PSD and 
likewise for the corresponding ACF. 


where we have introduced the root-mean-square arn- 
pl itude of our thermally driven motion as Xth = 
s/l/fimul = \JknT/rriuify In order to satisfy this nor¬ 
malization condition, we must integrate over the PSD 
once we have determined its functional form by evaluat¬ 
ing the convolutions found in Eq. (471. 
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Before we move on to the quantum PSD, we provide a 
brief remark in regards to carrying out the above proce¬ 
dure. As can be seen above, the calculations performed 
using this method become increasingly tedious as n be¬ 
comes larger, mainly due to the increasing number of 
convolutions. However, this complexity can be alleviated 
slightly by breaking up the convolutions into smaller cal¬ 
culations, allowing us to calculate our PSDs in an iter¬ 
ative manner which utilizes previous calculations. For 
instance, if we have already determined the second-order 
PSD, for which we need y(w) *%(w), we can convolve this 
quantity with x(w), or with itself, and use the result to 
determine the third- and fourth-order PSDs, respectively, 
reducing the number of convolutions needed. 


B. Quantum 


We now move onto calculation of the quantum PSD for 
x k ( t ). In this case, it is easier to focus on calculating the 
ACF, which can then be Fourier transformed to produce 
the corresponding PSD. Using Wick’s Theorem E3dB], 
we are able to determine the /cth-order ACF to be (see 

0 

R x k x k{t) = (x k (t)x k ( 0)) 

N 

= E Ac {x(t)x(0)) k ~ 2c (x(0)x(0 )) c , 

c—0 

(50) 


where 


A c 


(fc !) 2 

2 2c (c!) 2 (k — 2c)! ’ 


(51) 


and 


N = 


k/2 

{k — l)/2 


for even k, 
for odd k. 


(52) 


This quantity can be interpreted as the average number 
of quanta obeying Bose-Einstein statistics at a tempera¬ 
ture determined by /3. Combining the results of Eq. (53) 
with Eq. (501 we obtain 


N 


R x k x k ( t) — x 2k f A c (2 (n) + 1) 


c=0 


(55) 


x [((n) + l)e 


— iojQt 


(n) e 


,iu 0 t' 


k—2c 


Using the binomial theorem, we can instead write our 
/cth-order ACF in the form 


N k—2c 

Rx k x k (t) * Z pf£ ^ ^ B c d,C 

c =0 d —0 

where we have the new coefficient 


i(2c-\-2d— k)ujQt 


Bed, = 


( k\f (2 (n) + l) 2c ((n) + i) fc-2c-fl 


(n)° 


2 2 c (c!) 2 d! (fc — 2c — d)\ 


(56) 


(57) 


In this form, we can easily Fourier transform Eq. (56) to 
obtain the /cth-order quantum PSD 


iS fc rQ fc 


fc /j. k 


(t)e 

N k—2c 


iult dt 


M = 

= 27ra; 2 pf ^ Bcd6 


(2c 


c=0 d— 0 


2d — fc)wo), 
(58) 


where we have used the definition of the Dirac delta func¬ 
tion given by Eq. (All). This provides an expression for 
the /cth-order PSD for the undamped quantum harmonic 
oscillator for any positive integer k. 

In any realistic system, however, a non-zero amount of 
damping will occur as the oscillator radiates energy to its 
environment. To determine the the fcth-order PSD with 
damping included, we could in principle use Eq. (43) to 
calculate our PSD according to 


5^u k ij| k 


^=1 


(4*V) 4*V» ( 59 ) 


Therefore, by using Wick’s theorem, we have reduced 
the complex problem of finding a 2k term correlation 
function to evaluating the two term correlation functions 
under the sum in Eq. (50), which we determine to be (see 

H 


<*(t)*(0)> = [(<«> + 1) + <n> e***] , 

(x(t)x{t)) = (x(0)x(0)) = xl vf [2 (n) + 1]. 


Here we point out that in the first line we have obtained 
the expression for the ACF in the k = 1 (linear) case. 
In the above equations, we have introduced the thermal 
average of n for the harmonic oscillator, which is given 

by 


(n) 


1 

g/3?Kx>0 _ 1 ’ 


(54) 


where x^ioj) = T{x k (t)} = x(ui) * ... * i(w)/( 27 r) fc_1 . 
Here we have included a tilde over this PSD symbol to 
indicate that it is a quantum PSD with damping in¬ 
cluded. As we can see here, because we must work in 
the frequency domain for the input-output formalism of 
the damped harmonic oscillator, to determine the quanti¬ 
ties x^(u) we must compute k — 1 convolution integrals. 
Just like in the classical case, this leads to an increasingly 
complex problem as we increase k. 

Fortunately, by using a definition of the delta function, 
we have an alternate method by which we can include 
damping into the /cth-order PSD. In the case of small F, 
we can approximate the delta functions in Eq. (58) using 
Eq. (A12) to obtain 


4M 


1 

27r 


k r 


ur 


(fcr/2) 


2 ' 


(60) 














The subscript k is added here to differentiate between 
delta functions of different orders, as the half-width of 
the peaks of the PSD increases as kT/2. This effect is 
discussed in detail in Section [VI A| Using the expression 
in Eq. (60), we can write our fcth-order damped PSD as 


N k—2c 
2 k 


kTB cd 


s x k x k (to) = x 2k { y y - 

p ^ t'o (w + (2c + 2 d- k)oj o y + (kT/2Y 

(di) 

where we have a sum of Lorentzians instead of delta func¬ 
tions, effectively introducing damping into our quantum 
PSD. We will show below that for the k = 1 case, this 
result is exactly what would be obtained if we had in¬ 
stead decided to use the input-output formalism to in¬ 
clude damping in our system, justifying this simpler ap¬ 
proach. 

In concluding this section, we would like to point out 
that it is possible to use our result for the fcth-order ACF 
to determine the thermal average of x 2k (t). Taking t = 0 
in Eq. (501 we have 


~2 k 


>=*zp f (2< 


*> + l) 


k ( 2 k)\ / (H) \ ( 2 k)\ 


2 k k\ 


2 k kl 


(62) 


where ( H) = Hu>o{{n) + 1/2) is the average energy of the 
harmonic oscillator. This equation, unlike the classical 
analog, is valid for all temperatures as demonstrated by 
the fact that Eq. (49) is recovered by taking the high 
temperature limit k B T hu>o, for which 2 (n) + 1 ~ 
2k B T/Hu>o. 


The T = 0 limit of Eq. (62) can also be taken. Upon 
inspection of Eq. (54), we see that as T —► 0, (n) —► 0 
indicating that the oscillator is in its ground state, giving 


~2 k 


) = 


( 2 * 0 ! 

l zp{ 2 fcfc! • 


2k 


(63) 


We point out that this equation provides a quantum ana¬ 
log to Eq. (49), where we have taken x th x Z pf, as our 


system is purely driven by quantum fluctuations in the 
ground state as opposed to the classical thermal drive. 


V. FIRST-ORDER POWER SPECTRAL 
DENSITY 

Now that the framework for determining PSDs and 
ACFs for the harmonic oscillator has been laid out, we 
show that for k = 1 our formalism reproduces the well- 
known results of the first-order PSD for the position of 
the harmonic oscillator in both the classical and quantum 
regimes. 


A. Classical 

The linear displacement PSD for the classical damped 
harmonic oscillator is determined by taking k = 1 in 


Eq. (47), producing 


Sxx(m) — S FF (u}) |x(w)| — - 


Cth 

°FF 


m 2 [(w 2 - w 2 ) 2 + uj 2 T 2 } ’ 
(64) 

where we have used the generalized mechanical suscepti¬ 
bility found in Eq. (12). In order to determine the con¬ 
stant Sfip, we integrate S xx (oj) over all frequencies (see 
0 and use Eq. (36) to obtain 


Cth 

(x 2 >= * FF 


dui 


Cth 

°FF 


27tto 2 (w 2 — Wq) 2 + w 2 r 2 2m 2 WgP 


Inputting k = 1 into Eq. (49) we also have 


<x 2 ) = 


k B T 




This result, which can be written in the form 


<z 2 > = (V) = X k B T , 


(65) 


( 66 ) 


(67) 


is simply the equipartition theorem for the classical har¬ 
monic oscillator in thermal equilibrium at a temperature 
T, for which the average potential energy V is equal to 

k B T/2 |29| . 

By equating Eqs. (65) and (66), we find 
2mTk B T, which allows us to write the displacement PSD 
for the classical damped harmonic oscillator as 


Cth 

o F f 


Sxx(s*j') — 


2 Yk B T 


m [(to 2 — w 2 ) 2 + w 2 r 2 ] ’ 


( 68 ) 


This result agrees with that found in the literature US¬ 
ED, provided we incorporate the factor of 2 required 
when transferring between one- and two-sided PSDs. The 
above result could have also been obtained in a more 
straightforward route using the classical fluctuation- 
dissipation theorem GS3. which states that 


S X x(u) = —— Im{x(w)}. 


UJ 


(69) 


With the functional form of our PSD, we are now able 
to investigate some of its properties. First, since the 
signal is peaked at the resonance frequency and the PSD 
is an even function, we know that peaks exist at ui = ±wo, 
which results in 


5“ ax = S xx (±u 0 ) = 


2 k B T 


mu)Z 


(70) 


Another interesting parameter of the PSD is the width 
of the peak, which is closely related to the damping of the 
oscillator. Here, we consider the full width at half maxi¬ 
mum (FWHM) Aw. To determine this quantity we look 
for the frequencies Wi/ 2 at which S xx (u i/ 2 ) = <S£lf x /2, 
which leads to the quartic equation 


0 jj /2 + (r 2 - 2 wq) w\ /2 + w 2 (w 2 - 2r 2 ) = 0. (71) 
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Using the quadratic formula, the solutions to this equa¬ 
tion are found to be 

w i/2 = ±\Jul- Y -± yy'I 2 + ~ ±u} o ± (72) 

where we have made the high-Q approximation. The four 
solutions in the above equation correspond to two points 
on the sides of the two peaks at ±wo, which leads to a 
FWHM of Aw = T in the high-Q limit. 


B. Quantum 


Moving to the quantum regime, we now look to de¬ 
termine the first-order ACF and PSD for the position 
operator of the quantum harmonic oscillator. We have 


already calculated the first-order ACF in Section IV B 
where it was found to be (see the first line of Eq. (531) 


Rxx{t) = {x(t)x( 0)) 


= xtp { [((n) + l)e- i “ ot + (n)e i “° t }. 


(73) 


By Fourier transforming this ACF, or equivalently taking 
k = 1 in Eq. (581, we obtain 


S xx (u) = 2nxl f [((n) + l)<5(w - ui 0 ) + (n) 5(u + w 0 )]. 

(74) 

This produces the well-known expression for the first- 
order PSD for the position operator of the quantum 
harmonic oscillator [22]. This result is also verified 
by an independent determination using the fluctuation- 
dissipation theorem (see U>- 



0 


FIG. 1: The first-order damped quantum PSD, with each 
trace normalized such that its maximum at w = Wo is 1, vs 
frequency in terms of the resonant frequency. We have chosen 
a relatively low quality factor of 4 for this PSD, so that the 
amplitude effects are not obscured by the narrowness of the 
peaks. The different colors represent different average quanta 
with (n) = 0, 0.2, 0.5, 1, 2, 5, 100, as we move up in color 
from black to red. While the relative height of the peak at 
u) = wo is unchanged, the peak at ui = —u>o decreases with 
average number of quanta, demonstrating the asymmetry of 
the quantum PSD. In the extreme case of (n) = 0, the ui = 
—wo peak disappears altogether. 

The first-order quantum PSD is not symmetric about 
zero frequency as was the case for the first-order classical 


PSD. This asymmetry is visualized in Fig. [l] Physically, 
these two peaks correspond to two different processes. 
The negative frequency peak (w = — Wo) is associated 
with the annihilation/emission of a single quantum with 
frequency uiq. Alternately, the positive frequency peak 
corresponds to the creation/absorption of a quantum at 
Wo- In the context of optomechanics, these processes 
are strongly tied to Stokes/anti-Stokes Raman scattering 
whereby phonons can be created/annihilated via interac¬ 
tion with cavity photons [2]. Furthermore, the asymme¬ 
try of these peaks leads to distinctly non-classical effects 
at low phonon number, such as motional sideband asym¬ 
metry, which has recently been observed experimentally 

1231 EQ. 

It is also interesting to investigate the T = 0 limit of 
the above quantum PSD. This limit corresponds to the 
PSD of a quantum harmonic oscillator that is purely in 
its ground state, its motion arising solely from zero point 
fluctuations due to quantum noise. Taking (n) = 0 in 
Eq. (74), the quantum PSD becomes 


SxxM = 27ra z P f% - w 0 ). 


(75) 


In this limit, we completely lose the peak at w = — ojq 
due to the fact that in the ground state no quanta exist 
to annihilate. 

The above discussion on the physical significance of 
the quantum PSD was for the ideal case of zero damp¬ 
ing, leading to perfectly narrow peaks corresponding to 
quanta at two distinct resonance frequencies, ui = ±ojq. 
In a realistic system, however, damping will emerge, 
broadening these peaks and allowing for small deviations 
from this resonance frequency. We now look to include 
damping into our system by using the input-output for¬ 
malism outlined in Section [II C| In this case, we calculate 
the PSD using Eq. (43) to obtain 


Rxxi^') 




(x 1 (ui)x 1 (d)) did. 


(76) 


Utilizing Eq. (29), we find this damped PSD to be (see 

0 


Sxx(d) — rx zpf 


(n) + 1 


i) 


(ui — wq) + (r/2) (oj + uiq) 


( r / 2 ) 

(77) 

which also agrees to what is found in the literature [221 
123] , A T = 0 PSD corresponding to the ground state of 
a damped harmonic oscillator can also be determined by 
setting (n) = 0 in the above equation to obtain 


s°xx (w) = 


Ft 2 

1 •''zpf 


(ui — w>o) 2 + (r/ 2 ) 2 " 


(78) 


Both of these above results could have also been obtained 


by simply taking k = 1 in Eq. (611, justifying the method 


by which we obtained this expression. 

Finally, we can also find the thermal average (x 2 ) for 
the quantum harmonic oscillator using Eq. (42), where 
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we can integrate over either S xx (ui ) or S xx (uj) to obtain 

(H) 


(x 2 ) — x 2 p{ (2 {n) + 1) — 


raw n 


(79) 


consistent with k = 1 in Eq. (62). In the ground state, we 
then have ( x 2 ) = x 2 f such that the average value of the 
squared motion is the zero point fluctuation amplitude 
squared, as would be expected. This above equation can 
also be recast into 


1 2 
-mw 0 




(80) 


This can be interpreted as a sort of “generalized” equipar- 
tition theorem in which the average potential energy can 
be related to the average value of position squared, re¬ 
gardless of whether the drive results from thermal or 
quantum noise. As such, by taking the high temperature 
limit of this equation, the classical equipartition partition 
theorem given in Eq. (671 is recovered. 


C. Classical Correspondence 


For any quantum mechanical model, the correspon¬ 
dence principle tells us that the quantum result will re¬ 
produce its classical analog when the appropriate limits 
are taken. We will now show that for the first-order PSD 
for the damped quantum harmonic oscillator calculated 
above, we are able to retrieve the classical linear PSD in 
the limits of high temperature (IgT> fiwo) and quality 
factor (Q 1/2). This second condition must be taken 
as we have implicitly made assumptions of high-Q when 
introducing damping into the PSDs for the quantum har¬ 
monic oscillator. 

To begin, we see that in the classical limit we can use 
Eq. (54) to make the following approximation 


i) + 1 « (ra) 


1 


1 


_ = _ = k B T 

1 + /3fiwo — 1 (3huio hujo 


(81) 


Physically, this equation tells us that at high tempera¬ 
tures, the thermal energy of the resonator is broken into 
a large number of n quanta, each with energy Hojq, such 
that the ground state energy can be neglected. Remem¬ 
bering that x 2 p{ = 2 mu 0 ’ we can rewr it e Eq. |77l as 


£>XX (k-0 


rk B T r i i 

2mw„ _(w — w 0 ) 2 + (T/2) 2 (w + w 0 ) 2 + (T/2) 2 


2Tk B T 


m [(w 2 — uJq ) 2 + w 2 T 2 ] ' 


(82) 


The details of how the approximations were made to 
achieve the classical result are outlined in [Jj In the last 


line, we have retrieved the classical result of Eq. (681 that 
was determined in SectionlV A| Therefore, we have shown 


that the first-order quantum PSD calculated here satis¬ 
fies the correspondence principle in the region of interest 
surrounding the peaks at w = ±w 0 . This is illustrated in 
Fig. [2] for (n) = 1000 and Q = 1000. 



FIG. 2: Plot of the classical PSD (solid green) given by 
Eq. ([681 and the damped quantum PSD (dashed black) given 
in Eq. o in the regions surrounding the peaks at u = ±uio, 
exhibiting classical correspondence for the first order case. We 
have used (n) = 1000 and Q = 1000 such that the plots shown 
span 10r on either side of each peak. The PSDs are normal¬ 
ized such that their maximum values are 1 and the frequency 
is given in terms of the oscillator’s resonance frequency. 


VI. SECOND-ORDER POWER SPECTRAL 
DENSITY 

The linear PSD calculated in the above section can be 
used for situations in which the displacement of an oscil¬ 
lator is measured directly. However, there are situations 
where it is useful to measure the square of the position 
directly [3 IMS], in which case we need to consider the 
second-order PSD for the oscillator. In this section, we 
shall determine this quadratic PSD for both the quantum 
and classical cases. 


A. Classical 


We begin by calculating the classical second-order 
PSD, proceeding as we did in the previous section where 
we now use Eq. (47) with k = 2. The PSD in this case 
will be given by 


s x 2 * 2 (w) = S % F 2 IxM * y(w)| 2 , (83) 

the functional form of which is determined by the convo¬ 
lution 


xM * xM 


x(w')x(w 


uj') du'. 


(84) 


This integral can be computed using contour integration 
(see E to obtain 


xM * xM 


^47T* 

to 2 (w + IT) (w 2 — 4wq + 2 iujT) 


(85) 
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Inputting this expression into Eq. (83), we determine the 
unnormalized second-order PSD as 


first-order case. Evaluating the PSD at these two reso¬ 
nant frequencies, we find 


S x 2 x 2 


(w) = 


lQn 2 S t p-2 F 2 


to 4 (a , 2 + P 2 ) ^(w 2 — 4wq) 2 + 4w 2 r 2 ) 


■ ( 86 ) 


Upon investigation of this function, we see that it is 
peaked at ui = ±2a,o, as well as u> = 0. This is what we 
expect for the PSD of the squared displacement mi, 
the physical meaning of which will become more appar¬ 


ent when we look at the quantum case in Section VIB 


In order to properly normalize this second-order PSD, 
we must determine the value of S^ 2 p2 such that 


(x 4 ) =-!- / S x 2 x2 (uj)duj = 3 , (87) 

2W-oo V TOW o / 

where we have simply taken k = 2 in Eq. ( |49[ ). This 
is done by explicitly performing the integral by using 
contour integration (see 0 resulting in 


q2uq 

°x 2 x 2 


= S x 2 x 2(±2lo 0 ) 


(4a , 2 + 3T 2 ) (k B T) 2 
Tm 2 a, 4 (4cag + T 2 ) 


1 (k B T\ 
T \mul) 


( 93 ) 

where we have again made the high-Q approximation for 
the last step. Comparing these three maxima (see Fig.|3j) , 
we find that there is a global maximum at a, = 0 , with 
two local maxima at u> = A2a,o- As well, we find that for 
high-Q, the maxima are related by 5^2 = 4 S 2 2 ° 2 ■ 



(x*) = 


3TT 2 Sf 2F 2 


Tcj 2 m 4 (4a , 2 + 3T 2 )' 


( 88 ) 


Combining this with Eq. (87), we can solve for S F % p2 for 
which we find 


cth 

Op2 P 2 


Tm 2 (4a ; 2 + 3T 2 ) (k B T) 2 


(89) 


Putting this result into Eq. ( 86 ) we get the final form for 
the second-order PSD given by 


16T (4wq + 3T 2 ) {k B T) 2 

x2x2(UJ) m 2 . 2 (a , 2 A P 2 ) ((co 2 - 4a , 2 ) 2 + 4a, 2 T 2 ) ’ 

. (90) 

This equation can be simplified if we consider the high-Q 
limit, in which case we can express our PSD as 


*5 X 2 X 2 (^0 


64T (k B T) 2 

m 2 (a , 2 + P 2 ) ^(a ; 2 — 4a,g ) 2 + 4a, 2 T 2 ^ 


(91) 


FIG. 3: The second order classical PSD (solid green) in 
Eq. (|90|) and its quantum counterpart (dashed black) given in 
Eq. ( | 100[ ) in the regions surrounding the peaks at u> = 0, ±2o,o- 
Here we have used (n) = 1000 and Q = 1000 such that the 
plots are presented in a frequency window of 10P on each side 
of the peak as before. These PSDs are normalized such that 
their DC peak value is 1. 


We can also use these maxima to determine the 
FWHM at each of the peaks. Here we use a different 
approach than we did for the first-order case, due to the 
differing peak heights, as well as the fact that the denom¬ 
inator depends on the frequency to the sixth power. To 
simplify our calculations, we assume the high-Q limit to 
begin with and expand about a small deviation from the 
resonance frequency Slo, which at the half maximum we 
assume to be on the order of T. 

We begin by recalculating the width of the peaks in the 
first-order case to verify the effectiveness of this method. 
To do this we evaluate the first-order PSD at u = loq + Suj 
and equ ate it to half of its maximum value, which is given 
in Eq. (70), resulting in S xx (luo + Suj) = k B T/muj 2 Y. 
Evaluating this expression gives 


As we did above, we will now investigate the maximum 
values of the second-order PSD. The second-order PSD 
has three peaks, corresponding to three local maxima. 
Beginning by evaluating the peak at a, = 0, we find 


5; 


DC 


= s x 


»( 0 ) 


(4a , 2 + 3T 2 ) (k B T) 2 
r?n 2 a,g 


4 ( k B T \ 2 
T [mu 2 ) ’ 


( 92 ) 

where in the last step we have taken the high-Q approx¬ 
imation. 

We now perform the same calculation for a, = ±2a,o- 
Fortunately, due to the symmetry of the second-order 
classical PSD, both of these peaks will have the same 
maximum value, just as the a, = ±o,q peaks did in the 


(o,o A do,) — o,q A (ujq A $o,) T 2 — 2T“ 
S 4a,g(do,) 2 — CtlgT 2 = 0, 


(94) 


=k <5a, = A —, 


where going from the first step to the second step we have 
neglected all terms higher than second-order in Suj ~ F. 
From this result, we find the FWHM as Ao, = 2|do,| = T 
as we found above, verifying that in the high-Q case, 
the two methods give the same result. By using this ap¬ 
proach, we have implicitly assumed that we are only deal¬ 
ing with the a, = a,o peak, while ignoring the a, = —a,o 
peak. This is due to the fact that in this analysis, we are 
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only concerned with frequencies resulting from a small 
expansion around the peak of interest. However, this is 
inconsequential as the symmetry of the PSD ensures the 
negative frequency peak will have the same result. 

Now that we have verified the efficacy of this method, 
we apply it to the second-order PSD at w = 0. Evaluating 
our PSD at Soj and equating it to the half maximum at 
the DC peak we find 









[(<5o;) 2 + T 2 ] [((M 2 - 4wg) + 4(<5w) 2 r 2 

« 16wq((5u;) 2 - 16a^r 2 = 0, 

=> Soj = ±T. 


32r 2 a4 


(95) 


From this result, we find the FWHM at the DC peak to 
be Aojdc = 2|(5a;| = 2r, which is twice the value of the 
peaks in the first-order case. 

We now perform the same analysis for the peak at 
oj = 2loq (as usual symmetry ensures the same results at 
a i = —2u>o). Here we now expand about 2oco, such that 
we evaluate our PSD’s half maximum at oj = 2ojq + Soj, 
resulting in the expression 


[(2w 0 + M 2 + r 2 ] [ ((2wo + M 2 - 4wg) 5 

+ 4(2u; 0 + 5ui) 2 T 2 - 32r 2 a;^ (4wg + T 2 ) 

« 64wq(<5o;) 2 - 64a;^r 2 = 0, 

=> Soj = ±r. 


(96) 


Therefore, the FWHM of the peaks at w = ±2u;o is 
Ao; 2 W0 = 2\Soj\ = 2r, which is identical to what we found 
for the DC peak. 

It is interesting to note that the second-order PSD 
has peaks with twice the width of the first-order PSD. 
We will briefly investigate what this means physically. 
For the first-order PSD, the FWHM gives a measure of 
the spread of frequencies around resonance that a single 
quantum will have. The larger the width, the larger this 
spread. Now for the second-order PSD, we are looking 
at processes involving two quanta. Therefore, the spread 
of accessible frequencies doubles, as we combine the fre¬ 
quency distributions of each. This becomes mathemati¬ 
cally apparent by investigating Eq. (83), as the overlap of 


the mechanical susceptibility with itself in the convolu¬ 
tion integral reaches its half-maximum value when each 
first-order peak is ~ T away from the other, producing a 
second-order function with peaks of twice the width. See 
Fig. a for a more detailed explanation of this effect. 


B. Quantum 

We now look to calculate the second-order PSD in the 
quantum case. To do this, we follow a methodology sim¬ 
ilar to what we used to calculate the first-order quantum 
PSD by first finding the second-order ACF and using it to 


FIG. 4: A demonstration of how a convolution of two peaked 
functions with a width of F produces a peak with a width 
of 2r. In panels a)-e) the green curve is swept across the 
stationary blue curve while the red curve in f)-j) maps out 
the shared overlap area beneath the two functions, which is 
how we can conceptually think of convolution. In panel a) 
there is little overlap between the blue and green curve, so 
our curve in f) does not yet exist. In b), the peak of the 
green function has approached to within F of the blue. In 
this case, the overlap causes half of the green or blue curve 
to be integrated, such that in g) the red curve is at its half¬ 
maximum value. In c), when the blue and green functions are 
perfectly overlapping, the red curve in h) is at its maximum 
value. In d) we have the opposite case to b), where the green 
curve is now receding and the red curve in i) is decreasing 
through its half-maximum. Finally in e), the green curve has 
left the frame and the majority of the overlap with the blue 
function has occurred. The end result is the red function in 
j) with twice the width of the two peaks that were convolved 
to produce it. As a final note, we mention that this process 
can be performed any number of times, convolving k curves of 
width r a total of k — 1 times to produce a final distribution 
with width kT. Thus we have justified our choice of half-width 
kT /2 in Eq. (601. 


find the corresponding PSD. Starting by determining the 
second-order ACF, we take k = 2 in Eq. (551 to obtain 


- T 4 

— "‘'zpf 


Rx*x*(t) = (x 2 (t)x 2 { 0)) 

= 4>f ( 2 [«"> + !) + <«> e***] 2 + [2 (n) 2 + 1 

2 ((n) + l) 2 e~ 2iUot + 2 (n) 2 e 2i “°* 

+ 8 (n) ((n) + 1 ) + 1 . 

(97) 

By Fourier transforming the above equation we obtain 
the second-order PSD of the quantum harmonic oscillator 


S x 2 x 2 {oj) = 2nXz p f 2 ((n) + l) 2 S(oj - 2u ; 0 ) 

+ 2 ( n ) S(oj + 2ojq) + (8 ( n) ((n) 1) + 1) (5(o;) 


(98) 


which could also have been obtained by taking k = 2 in 
Eq. (58). This result agrees with that found in Eq. (S3) of 
El and is also verified using the fluctuation-dissipation 
theorem (see [K| ) . 

Eq. (98) also exhibits asymmetric qualities, similar to 
the first-order case, which are displayed in Fig. [5] As 
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FIG. 5: The second-order quantum PSD, with each trace nor¬ 
malized such that the peak at to = 2too is 1, vs frequency in 
terms of the resonant frequency. The color scheme is the same 
as in Fig. [l]so that (n) = 0, 0.2, 0.5, 1, 2, 5, 100, as we move 
up in colors and a quality factor of 5 is chosen for clarity. In 
this case, the relative height of the peak at to = 2too stays con¬ 
stant, while the other two peaks decrease as (n) is reduced. 
Notice that in this case, when (n) = 0, the to = — 2wo peak 
vanishes, while the DC peak still has a finite height, albeit 
significantly reduced. 


we would expect from the classical case, this function 
contains peaks at u = 0, ±2wo. Furthermore, by in¬ 
specting the coefficients of (n) 2 in each term, we find 
SB£ a = 4S 2w °2 for the large (n) case associated with the 

X Jb XX ^ ' ' 

classical regime, agreeing with what we found in Section 

mb 

Assigning a physical interpretation to the above peaks, 
two quanta are annihilated at the to = — 2too peak, while 
we have creation of two quanta at the to = 2wo peak. 
However, we now have a new DC term at ui = 0. Phys¬ 
ically, this peak corresponds to the second-order process 
by which either a quantum is created then annihilated, or 
annihilated then created, with no net change to the sys¬ 
tem. These processes are unique to a nonlinear system, as 
we require a two-step procedure, which is prohibited for a 
linear system. In cavity optomechanics, this term causes 
a DC shift in the optical cavity’s resonance frequency in 
proportion to the number of phonons in the mechanical 
resonator, providing an avenue by which we can perform 
a QND measurement of these quanta j6\ In addition 
to this effect, cavity optomechanical systems also exhibit 
other second-order effects corresponding to two-phonon 
processes, leading to phenomena such as mechanical cool¬ 
ing/squeezing [5], as well as optomechanically induced 
transparency mm- 

As we did in the first-order case, it is interesting to 
investigate the T = 0 limit of the second-order quantum 
PSD. Taking (n) = 0 we are left with 

S° 2 x 2 (uj) = 2TTxi pi [2S{u) - 2to 0 ) + <5(w)]. (99) 

Again, processes involving the initial annihilation of 
quanta vanish at T = 0, leaving the two remaining peaks 
at w = 2ujq and to = 0, corresponding to the creation of 
two quanta and the creation of a single quantum followed 
immediately by its annihilation. What is surprising is 


that the relative height of the DC peak and the 2too peak 
has decreased eight-fold, as the DC peak is now half of 
that at 2ujq. 

We can also determine the second-order PSD of the 


damped harmonic oscillator by taking k = 2 in Eq. (61) 
for which the result is 


S x 2 x 2(to) = 2Txt p{ 


2((n) + l) 2 


2 (n) 2 


(to — 2wo) + r 2 (to + 2 wo) 
8 (n) ((n) + 1) + 1 


+ r 2 


r 2 


( 100 ) 


The zero temperature limit of the above equation can be 
obtained by taking (n) = 0 as usual, resulting in 


^x 2 x 2 ( w ) — 2FXz pf 


1 


(to — 2 wq) + r 2 to 2 + 


r 2 


( 101 ) 

It is also interesting to investigate the thermal average 
(£’ 4 ) for all T. Inputting either the undamped PSD of 
Eq. ([981) or the damped PSD of Eq. (|100|) into Eq. (|42|) 


as we did in the first-order case we obtain 


(£ 4 ) = 3x 4 pf (2 (n) + l) 2 , 


( 102 ) 


consistent with the k = 2 case of Eq. (62). For T = 0, this 
reduces to (x A ) = 3cc 4 pf , which is in direct correspondence 


to the classical case given in Eq. (87), where we have 


replaced the thermal drive with a quantum one. 


C. Classical Correspondence 


We conclude this section by ensuring that our second- 
order quantum PSD obeys the correspondence principle. 
Beginning with Eq. (100) and using the same approxi¬ 
mations as we did for the first-order PSD. we find (see 

0 


- r ( k BT\ 2 

^x 2 x 2 (^) ~ r i 2 ) 

1 , 1 4 
x [(to - 2w 0 ) 2 + r 2 + (to + 2to 0 ) 2 + r 2 + to 2 + r 2 _ 

64F (k B T) 2 

m 2 (to 2 + F 2 ) ((w 2 — 4wg) 2 + 4w 2 r 2 ) 

(103) 


which matches the expression we found in Eq. © such 
that classical correspondence is again satisfied. This is 
presented visually in Fig. [3j 


VII. CONCLUSION 

We have presented a method to calculate a general 
PSD for the classical and quantum harmonic oscilla¬ 
tor, corresponding to any power of its position that is 
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a positive integer. We then investigated the experi¬ 
mentally relevant cases of k = 1 and k = 2 associ¬ 
ated with the linear and quadratic PSDs, respectively. 
The expressions for the first-order PSD are well-known 
[HH25] and are presented here to verify our general 
model. For the case of the second-order PSD, a number 
of expressions useful in the field of quantum measure¬ 
ment were calculated, allowing researchers to ascertain 
whether or not a device is suitable for QND measure¬ 
ments Both of these results were 

found to agree with an independent determination rely¬ 
ing on the fluctuation-dissipation theorem. 


Though this document focused largely on the first- 
and second-order cases, as these are the regimes that 
are closely linked with experiment, it is possible higher- 
order PSDs may become useful in the near future. For 
instance, quartic or fourth-order coupling has already 
been achieved |4] and is proposed as a method to gen¬ 
erate Schrodinger cat states [521155] . For this reason, the 
results presented in this paper provide a useful tool by 
which a theoretical nonlinear PSD can be determined and 
fit to experimental harmonic oscillator position spectra 
of any order. 


Appendix A: Definitions 
1. Fourier Transforms 


In this document, we choose to define our Fourier transform for an arbitrary, time-dependent signal a(t) as 

/ OO 

a(t)e lult dt = a(w), (Al) 

-OO 

with the inverse Fourier transform being defined as 

1 r°° 

J r_1 {a(w)} = — / a{uj)e~ lut du = a(t). (A2) 

2tt J- oo 

With these definitions, we find the following useful property of the complex conjugate of the Fourier transform 


/ /»oo \ * /*oo 

a*(uj) = ( / a(t)e lult dt) = / a(t)e~ lut dt = a(—uj). 

\J — oo / J —oo 


(A3) 


Here we have assumed that the time-varying signal a(t) is real. We also have the expression for the Fourier transform 
of the nth time derivative of a(t) given by 


F { } = H w )" F{a{t)} = (—iuj) n a(u). 


(A4) 


2. Convolution 


For an arbitrary variable u, the convolution of two functions g(u) and h(u) is defined as 

/ OO 

g{u')h(u — u') du ', 

-OO 


(AS) 


and represents the measure of the area shared by the two functions, as we translate one across the other. A useful 
property of this operation is that the Fourier transform of a convolution of two signals in the time domain is the 
product of their Fourier transforms in the frequency domain. That is to say 


T{g{t) * h(t)} = g(uj)h{uj), 


(A6) 


where g(u>) = T{g{t)} and h(cu) = J-{h(t)} are the Fourier transforms of the signals g(t) and h(t). 

We can also use the convolution integral to express the Fourier transform of a product of two functions in the time 
domain as the convolution of their frequency domain representations. For our Fourier transform definitions, this is 
given by 

F{g{t)h(t)} = * h{u). 


(A7) 
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Finally, this property can be extended to a product of k functions in the time domain giving 

F{9i{t)9i{t)-9k{t)} = ( 27 r )fc-i g1 ^ * 9 * SfcM> 
where the ellipsis (...) is used to indicate that there are k terms in the sequence. 


(A8) 


3. Dirac Delta Function 


Conventionally, the one-dimensional Dirac delta function is defined as 


S(u — v!) 


oo if u = u', 
0 if it 7 ^ u', 


such that 


J f{u)5{u — u) du 


/(«)> 


(A9) 


(A10) 


provided that the integral contains u = v! in its range of integration. 

This is not the only way to express the Dirac delta function, however, and here we introduce two alternate definitions 
which are used above. The first is given by 


5(u — u') = 


2ir 


0 iv{u-u') dv _ 


(All) 


This definition proves to be very useful when performing Fourier transform integrals throughout this paper. 
The second definition arises from taking the limiting case of a Lorentzian function, such that 


1 


S(u — u) = lim — 


e->o 7 r (it — u'Y + e 2 ’ 


(A12) 


where e is the half width at half maximum (HWHM) of the peak. This equation provides one of the simplest ways for 
introducing width to peaks that are infinitesimally narrow ( i.e. including damping where there was initially none). 


Appendix B: Contour Integration 


Contour integration is a powerful method of integration by which functions containing complex poles can be easily 
integrated. This technique hinges upon the Cauchy’s residue theorem, which states [34] 


N 

f g{z)dz = 2ni V] Res (g, zi), 

J c 1=1 


(Bl) 


where the integral is performed over a positively oriented (counter-clockwise) closed contour C and Zi are the N poles 
of g(z) enclosed by C. As well, we have introduced the residue of g(z) at Zi denoted by Res(g, Zi). For a function that 
can be expressed as 


g(z) 


0(z) ? 

z - Zi 1 


we can determine its residue at Zi as 


Res (g,zi) = 0{zi). 


(B2) 


(B3) 


Often, a convenient choice for the contour C is a semicircle with infinite radius that extends over the region of the 
complex plane with Im{z} > 0. In this way, the portion of the semicircle that runs along the real axis stretches from 





16 


—oo to oo. If we assume that \g(z)\ falls off faster than 4 as z —> oo, the portion of the integral performed on the 
curved part of the contour will be zero, allowing us to write 


£g(z) dz 


g(z) dz 


N 

2m Y R.es(g, zi ). 

i=i 


(B4) 


Our choice of contour is such that we are now only concerned with the poles in the top half of the complex plane. We 
could have alternatively chosen the semicircle enclosing the poles in the bottom half of the complex plane, provided 
we account for the minus signs that will arise due to our differing contour orientation. The integral given in Eq. (B4) 
is useful in computing a number of quantities for PSDs. 


Appendix C: Thermal Average of x 2k for the Classical Harmonic Oscillator 

In order to normalize the PSDs calculated in this document, it is important to know the thermal average of the 
position to even powers, that is (x 2fc ). For a one-dimensional, classical system, the thermal average of a quantity A 
is given by 


(A) = 


JI Aer f3H dx dp 
jIe~ pH dxdp 


(Cl) 


where /3 = l/k B T and H = H(x,p) is the Hamiltonian of the system as a function of position x and momentum p 
and the integrals are performed over the entire region of each corresponding phase space. 

In order to calculate the therm al av erage of x 2k (t) for the one-dimensional harmonic oscillator, we input the 
Hamiltonian from Eq. ^ into Eq. (Cl) with A = x 2k resulting in 


<* 2fc > = 


/ OO /»00 2 /*00 

x 2k e~ipm U 2 0 x 2 dx e -^r dp / x 2k e~^ mu, ° x2 dx 

-oo J —OO J —oo 


dx 


tip 2 , 

e 2 ™ dp 


0 -\pmuilx 2 dx 


(C2) 


To evaluate the second line of Eq. (C2), we calculate the integral in the numerator, which can be shown (by induction) 
to be 


X ih e -\f>mu, 2 0 x 2 dx = 


(2 k)\ 


( 2f3mujQ) k k\ 7_ c 


e -2t>mu> 0 x dx _ 


(C3) 


We can therefore input Eq. (C3) into Eq. (C2) to obtain 


/ 2k\ ( 2fc ) ! _ 2fc ( 2fc ) ! 

\ ' (2/3 muj 2 ) k k\ th 2 k k\ 1 


(C4) 


where Xth = \A/ /3mujQ = \JksT/rruxi^ is the root-mean-square amplitude of our classical thermally driven motion. 
Inputting k = 1 and k = 2 we obtain 


and 


<x 2 > 


1 

/3towq 


k B T 

mu>Q 


L th> 


(C5) 


<x 4 > = 


(/Smwg ) 2 


= 3 


knT 


muif 


= 3x 


th’ 


(C6) 


which are of special interest for this document as they are required to properly normalize the first- and second-order 
classical PSDs. 
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Appendix D: Wick’s Theorem 


Wick’s theorem [271 GB] is a powerful operator identity that is often used in quantum field theory to simplify the 
products of creation and annihilation operators. Here we show how it can be used in the context of determining the 
fcth-order ACF for the position operator of the quantum harmonic oscillator. 

In order to properly introduce Wick’s theorem, we must first define a number of ordering operations which can be 
performed on products of quantum mechanical operators. The first such operation is known as time-ordering and 
will be enforced using the time-ordering operator T. When applied to a product of time-dependent operators, time¬ 
ordering ensures that operators with a larger time argument appear left of those with a smaller one. For instance, 
when applied to the product of bosonic operators x(ti)x(tj ), the time-ordering operator produces 


T{x(U)x(t.j)} 


x(ti)x{tj) for ti > tj, 
x(tj)x(ti) for ti < tj. 


(Dl) 


We can also introduce an anti-time-ordering operator 7”, which has the opposite effect of the time-ordering operator. 
That is to say the anti-time-ordered product is organized from left to right by ascending time arguments. Applying 
the anti-time-ordering operator to the example given in Eq. (|D1|) now instead gives 


T{x(tj)x(tj)} 


x(tj)x(ti) for ti > tj , 
x(ti)x{tj) for ti < tj. 


(D2) 


Next, we introduce the concept of normal ordering. This operation is performed by the normal ordering operator 
A7, which takes an arbitarily ordered product of creation and annihilation operators and arranges them so that all of 
the creation operators are on the left of the annihiliation operators. A simple example of this process is the normal 
ordering of the product bb^bb\ which is given by 

J\f{btfbP} = tftfbb. (D3) 


This operator ordering has the very useful property that when acting on the ground state, it always produces zero as 
a result of the rightmost annihilation operator. 

Finally, we introduce the contraction of two operators, which for the case of two position operators is given by 


x{ti)x(tj) = T{x(U)x(t.j)} - N{x[ti)x{tj)}. 


(D4) 


Using the expression for x[t) given in Eq. (24), we can explicitly calculate the contraction above as 


x{U)x{tj) = a4f x 


e b) for ti > tj , 

e -iui 0 (tj-ti) for ti < tj. 


(D5) 


where we have used the commutation relation [b. U] = 1. It is interesting to note that the contraction shown here is 
simply a complex number (not an operator) independent of the number of quanta in the system. 

We now have the machinery required to properly present Wick’s theorem, which states that a time-ordered product 
of operators comprised of creation and annihilation operators can be represented as the normal ordering of the product, 
plus a sum over all possible unique contractions of the operators within the product. For the position operators of 
the quantum harmonic oscillator, this is mathematically expressed as [28] 


T{i(tAr)i(tAr-i)...i(t 2 )^(U)} = JV{x(t N )x(t N -!)...x(t 2 )x(ti)} + ^ N {x(t N )x{t N -i) ...x(t 2 )x(ti)} 

single 

+ ^2 J^{x{t N )x{t N -i)...x(t 2 )x{ti)} + ... + y^ j JV{x(t N )x(t N - 1 )...x(t 2 )x(t 1 )}, 

double all 


(D6) 


where the sums are carried out by performing all of the specified unique contractions (single, double, triple, etc.) and 
N = 2k is the number of position operators we are considering. 

We now look at how this theorem can be used to reduce the correlation functions comprised of 2k operators (known 
as 2fc-point correlators) used to calculate the fcth-order ACFs to a sum of products of two-point correlators. We begin 
by introducing the following useful identity comprised of correlators of normal ordered position operators 

{AT{x{t N )x(t N - 1 )...x{h)x(t 1 )}) = ^(7U{i(t i Ji(t Jfc )})...(A7{i(t ll )x(t il )}), 

TO 


(D7) 
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where the indicates a sum over all unique products of time-ordered two-point correlators on the RHS of the 


TO 


equation, that is t* > tj for all sub-indices on i and j and there are no repeated terms. We also note that this identity 

/V /S A. * k 

hinges on the fact that we can write ((b^) k b k ) = k\ (b^b) . 

Next we rearrange Eq. (D4) and take the thermal average resulting in 


{T{x(ti)x(tj)}) = {.Af{x{ti)x(tj )}) + (x{ti)x(tj)). 
We can now imagine taking the product of k of these correlators which results in 


(D8) 


(x(t N )x(t N - 1 ))...(x(t 2 )x(ti))) = (Af{x(t N )x(t N -i)})...(J\f{x(t 2 )x(t 1 )}) + ^2{x(ti)x(ti- 1 )) J|(A/'{x(t;)^(^-l)}) 

i 


+ '^2 l (x{ti)x(ti-\))(x{t j )x(tj-\)) (Af{x(tl)x(ti- 1 )}} + ... + (x{t N )x(t N -i))...(x(t 2 )x{ti)), 

ij l&,3 


(D9) 


where we have taken t jv > tjv_i > ••• > t 2 > t\ such that we can drop the time-ordering operators on the LHS. In other 
words, the product on the LHS is equal to a product of k normal ordered two-point correlators, plus terms containing 
single contractions, plus terms containing double contractions, etc., all the way until we reach the term containing 
k contractions. Finally, we sum over all unique combinations of Eq. (D9) such that each two-point correlator in the 
product on the LHS is time-ordered to obtain 


^(x{ti k )x(t jk ))...(x(t ll )x(t h )) 


TO 


= ^(A/"{®(i ifc )®(tjj})...<A/'{^(ti 1 )i(t jl )}) +^2^2{x(t irn )x(t jrn )) 

l^m 


TO 


TO m 


+ J2J2^ ti ^^tj m )){x{t i Jx{tjJ) {N'{x(ti)x(t l -. l )}) + - + ^2(x(t ik )x{t jh ))...{x{t il )x(t jl )) 

TO m,n l^m,n TO 

= (Af{x(t N )x(t N - 1 )...x(t 2 )x(t 1 )}) + ^2 (J^{x(t N )x(t N - 1 )...x(t 2 )x{t 1 )}} + ^2 (■^f{x{t N )x{t N - 1 )...x(t 2 )x{t 1 )}) 


single 


double 


+ ...+£w x(fAr)£(tAr-l)...x(t 2 )£(tl)}), 


all 


(DIO) 


where we have used the identity in Eq. (D7) with the same definition of as above, as well as the fact that the 


TO 


contraction of two position operators is a complex number independent of the number of quanta in the system such 
that it can be pulled outside the correlator. By inspecting the final result of Eq. (DIO), we see that by summing over 
all unique combinations of k time-ordered two-point correlators, we obtain what we would get if we took the thermal 
average of Eq. (D 61 . Therefore, this allows us to use Wick’s theorem to write a 2fc-point position operator correlator 
as 




TO 


(Dll) 


It is this result that allows us to break our fcth-order ACFs into sums of products of two-point correlators. 

Before we go into more detail on how to calculate fcth-order ACFs using this result, we present a brief example to 
elucidate the above method by showing how it can be applied to reduce the four-point correlator (x(t 4 )x(t 3 )x(t 2 )x(ti)} 
to a sum of three products of time-ordered two-point correlators of the form (x(ti)x(tj)). For this example, we assume 
that t 4 > t 3 > t 2 > ti such that we do not need to explicitly write out the time-ordering operator. Using Wick’s 
theorem, we can write this four-point correlator as 


(x(t 4 )x(t 3 )x (t 2 )£(U)} = {M{x{t 4 )x(t 3 )x(t 2 )x(t-i)}) + (A /’{x(t 4 )x(t 3 )x(t 2 )x(t 1 )}) + (Af{x(t 4 )x(t 3 )x(t 2 )x(ti)}} 


+ (A/'{£(f 4 )£(t 3 )£(l 2 )i(tl)}) + (N {x(t 4 )x{t 3 )x{t 2 )x(ti)}) + (N {x{U)x(t 3 )x(t 2 )x(t 3 )}) + (Af{x(t 4 )x(t 3 )x(t 2 )x(ti)}) 

(D12) 


+ (W{£(f 4 )£(As)£(l 2 ):E(U)}) + {N {x{t 4 )x(t 3 )x(t 2 )x(ti)}) + (A/'{i(t 4 )i(t 3 )x(t 2 )i(ti)}). 
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Alternatively, we could have used Eq. (Dll) to write our four-point correlator as 


(x(t 4 )x(t 3 )i(t 2 )i(*i)) = {x(t4)x(t 3 )){x(t 2 )x(ti)) + {x(U)x(t 2 )){x(t 3 )x(ti)) + (^(t 4 )®(<i))(®(^3)^(ii)), (D13) 


where we have been careful to ensure that each of our two-point correlators are unique and time-ordered. Evaluating 
this correlator using either Eq. (D12) or Eq. (D13l, we obtain the identical result, namely 


(s(t 4 )s&(t 3 )x(t 2 )*(ii)) = xi pi 2 (n) 2 (e -i “ o( * 4+t3- * 2_tl) + e i “ o( * 4+t3_t2 "* l) + e -^o(t4-t 3 +t 2 - tl ) + 


| g — iuJo(t4 —t 3 — t 3 + ti) gicUo(t4—1 3 — ^2-t-q)^ ^^—i(^o(.t4+t 3 — t2 — tl) — C^4 — 1 3 +^2 ~ 1 1) 

g4LJo(*4 — t 3 +t2 — ti) 2 e ~ iul o{t4-t 3 -t2 + t 1 ) 2^ UJ o(t4 — t 3 — t2+tl)^j g — IUJq (t4 — t 3 +t-2 — 1 1) 2^>~^ UJ o(^4+t 3 —t2 — ti) 

(D14) 


We note that the above result is consistent with that found in Eq. (97) for the second order ACF, provided we take 
t -4 = t 3 — t and t 2 = ii = 0. 


Continuing on to calculate the fcth-order ACF, we now imagine replacing fjv = t^-i — ■■■ = tN-k+i = t and 
fjv-fc = tjv-fc-i = = ii = 0 in Eq. (Dll), which casts the correlator on the RHS into the form identical to that 

of Eq. (50). At first glance, it would seem that due to the implicit assumption of time-ordering in Eq. (Dll), this 
association is only valid for t > 0, where we are ensured that our fcth-order ACF is time-ordered. However, though 
we have been assuming time-ordering throughout this section, the exact same results can be shown to be true for 
anti-time-ordering (ze. replace all time-ordering operators with anti-time-ordering operators). Therefore, this result 
is also valid for t < 0. The trivial case of t = 0 can be handled by the continuity between these two domains. Our 
final result then becomes 


N 

(x k (t)x k ( 0)) = ^2 A c {x(t)x( 0)) fc_2c (z(0):e( 0)) c , 

c=0 


(D15) 


where A c is given by Eq. (51) and counts the terms of differing orders in the two-point correlators used here. 


Appendix E: Calculation of First-Order Quantum Autocorrelation Functions 


Here we calculate the two-point correlators comprised of x(t) and x(0) found in Eq. (50) which are used as the 
building blocks to determine any general ACF for the position of the harmonic oscillator to the fctli power. To calculate 
these three correlators, we use Eq. (24) to obtain 


<£(t)4(0)> = x 2 p{ ((66) e~ iu,ot + (66+) e -^ ot + (S+S) e iu>ot + <S t S t ) , 
(i(t)i(t)) = x 2 zp[ ((66) e~ 2iuiot + (66+) + (6+6) + (6+6+) , 

(z(O)x(O)) = x 2 zp{ ((66) + (66+) + (6+6) + (6+6+)) . 


(El) 
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To determine the correlators of the ladder operators found in the above equation, we use Eqs. (38) and (39) along 
with the properties of the ladder operators given in Eqs. and to find 


@b) = 


P e /3E " (n\ 6^6 \n) ^ e l3Erl n(n\n) y^ne 0En 


E' 


,-/3 E„ 


E< 


,~PE n 


E< 


~PE„ 


= {n) 


m = 


e (n| 66^ |n) E e ( n + 1) ( n \ n ) P(n + 1)< 


-0E„ 


E« 




E< 


,-/3S„ 


E< 




— (n) + 1 , 


P e /3E " (n| 66 \n) E e V n ( n — 1) (''T-p — 2) 


(E2) 


( 66 ) = 


E- 




E< 


~0E„ 


= 0, 


<StSt) = 


yy e (n| |n) E e l3En n ( n + 1) ( n l n + 2) 


E< 




E« 




= 0. 


As expected, only the correlators with one creation and one annihilation operator are nonzero. We have also introduced 
(n) which is the average thermal population of quanta and is given by 


(n) = 


E ne 

n —0 


—/3 hu) o ( n +1 / 2) 


E’ 

n —0 


0 —phujQn 


^ e -phu> 0 (n+l/2) y" e 


—(3 hujQn 


(Er?®o) 5 

l 

l_ e -tmoj 0 
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e /3hui 0 _ ^ ’ 


(E3) 


n—0 


n—0 


where in the above equation we have used Eq. to input an expression for along with the following sums 

OO 

E»* = V. 

f—' 1 — k; 


fc=o 


(E4) 


E*»‘ = V 


fc=0 


(1 — u >) 2 ’ 


It should be noted that the value we obtained for (n), known as the Bose-Einstein occupation factor, is exactly what 
we would expect for the average occupation number for a thermal distribution of bosons. 

Using the relations given in (E2), we can then express our first-order position correlators as 


(mm) = 4>f [«"> + !) e-“ ot + in) e***] , 

(mm) = (mm) = x 2 zp{ [2 (n) + 1 ]. 


(E5) 


These correlators are the ones given in Eq. (53), which are used to determine the fcth-order quantum ACF. 


Appendix F: Contour Integration of the First-Order PSD 


Here we use contour integration to calculate the integral given in Eq. (|65|. We begin by factoring the denominator 
of Eq. (64) in terms of its complex zeros using the quadratic equation, allowing us to write 




Cth 

o F f 


77V 


w - i (zr + v^-r 2 )] [ w - | (*r - y/4w* - r 2 ) 
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u +1 (iT - v/4^-r 2 )] [ w +1 (»r + - r 2 ) 
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We now have our linear displacement PSD in a form similar to that given in Eq. (B2), allowing us to easily compute 
the residues corresponding to the different poles. Inspecting the above equation, we can read off these poles as 
to = ±| (p ± \/4w q" — r 2 ). Since we are in the high-Q limit, we are assured that the square root quantities will 
be real, as 4wg > r 2 . Therefore, only two of these poles will be in the top half of the complex plane, namely 
ui = \ + ^/4ojg — and W 2 = \ — T 2 ) ■ We can then write the integral in Eq. (65) as 


1 r°° _ 

(x ^ 7^ J S xx {^) diO — i [Res(*S , ;c;r , CUi) T Res(*5 f ;ca :,CU2)] . 


(F2) 


Computing each of these residues separately using Eq. (B3), we find 

R-es(S' a;x , wi) 


Cth 

°FF 


Res(Sica,, W 2 ) 


im 2 r^Aco 2 - r 2 (*r + V4w 0 2 - r 2 ) ’ 

_ Qth 
°FF 


(F3) 


im 2 t^Alo 2 - r 2 (*r - v'^o 2 - r 2 ) 


Putting these results into Eq. (F2), we obtain 


<* 2 > = 


oth 

°FF 

2iti 2 u}qT 


(F4) 


Combining this result with Eq. (66), we can now determine Sftp, allowing us to properly normalize the first-order 
classical PSD. 


Appendix G: Contour Integration to Determine the Functional Form of the Second-Order PSD 


In order to determine the second-order displacement PSD, we need to perform the convolution integral x(w) *y(w). 
To evaluate this integral we follow a strategy similar to the previous section, in which we express the convolution 
integral given in Eq. (84) as 


/ oo ^ roo 

x(w')x(w - w') du}' = — / y(w, a/) du/, 

-00 m J- 00 


(Gl) 


where y(u>, u/) is a function expressing the integral in terms of its poles and is given by 
y(uj,uj') = T 


'w +1 (*r + V4wg-r 2 )] p +1 (*r - - r 2 ) 


(G2) 


x -F 


h (* r + \/4wg-r 2 )] p- w -|(iT-V4o; 2 -r 2 ) 


In this form, the poles are easily read off at <*/ = — | ^*P ± \/4wg — T 2 j and w' = u + | ^?T ± — T 2 ^. Again, 

choosing the semicircle that covers the top half of the complex plane, we concern ourselves with the two positive 
poles enclosed in this region, namely = u + | + y/4wg — T 2 ^ and lu 2 = w + \ ^*T — \J — P 2 ^ . We can then 


27ri 


express our convolution integral as 

xM * xM = ^ [Res(j/, w'J + Res(y,u4)] • 
m z 

We calculate these residues separately, obtaining 

Res(y,wi) = 


(G3) 
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Res(y, W 2 ) = 
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Inputting these two residues in Eq. (G3), we get 


xM * xM = 


—47 ri 


m 2 (tjj + ?r) (w 2 — 4 ujg + 2icoT) 


(G5) 


This is the expression that allows for determination of the functional form of the second-order PSD for the classical 
damped harmonic oscillator. 


Appendix H: Contour Integration of the Second-Order PSD 


We use contour integration one last time in order to determine the normalization constant S ^ 2F 2 for the classical 
second-order PSD. Following the same procedure as in the last two sections, we begin by writing the second-order 
PSD given in Eq. ( 86 ) in terms of its poles 


S x 2 x 2 (u) = 


167T 2 S^ 2i r 2 


(uj + *r) (w — iT) + *r + a/4wq — r 2 ) 


(HI) 


+ *r — yj 4wq — r 2 ) (u> — *r + \J 4wg — r 2 ^ (w — *r — \J 4wg — r 2 ^ 


Upon inspection of this equation, we can see that the poles are given by oj = ±*r and lo = iiT ± \J 4ajg + T 2 . We 
now look to compute the integral given in Eq. (871. Using our usual method, we concern ourselves with the three 
poles in the upper half of the complex plane, given by wi = *r, 0 J 2 = *T — \/ 4 cOq — T 2 and w 3 = *r + \/4uJq — F 2 . We 
can then express the integral in question as 


1 r°° 

(x 4 ) = — J S x 2 X 2 (uj)duj = i [Res(5' a; 2 a .2, oq) + Res(S' a; 2 a 2 , w 2 ) + Res^^a, w 3 )] . 
Computing each of the residues separately we find 
Res (5^2 3 , 2 , 
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(H3) 


Inputting these residues into Eq. (H2), we obtain 

<* 4 > = 


8K 2 S% F 2 

Twoto 4 (4^0 + 3r 2 ) ’ 


(H4) 


from which we can determine S mm using Eq. (871, allowing us to normalize the second-order classical PSD. 


Appendix I: Calculation of the First-Order Damped Quantum PSD using Input-Output Theory 

In this section, we use the expressions for a harmonic oscillator coupled to an external heat bath given in Section 


IIC to calculate the first-order PSD for a damped harmonic oscillator. Inputting Eq. (291 into Eq. (76), we can write 

dui 1 


our damped PSD as 
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where we have taken (& 7 (w)& 7 (w')) = (bj) (w)M) (w'H = 0 due to the last line of Eq. (30). Using the other two relations 
of Eq. (301, we can perform the integrals in Eq. © to find 


S X x(u) — Ta; Z pf 


(n) + 1 


+ 


i) 


(w - w 0 ) 2 + (r/2) 2 (w+w 0 ) 2 + (r/2) 2 


( 12 ) 


Here we have taken nb(wo) = (n) due to the fact that in thermal equilibrium the oscillator and heat bath will be at 
the same temperature. 


Appendix J: Classical Correspondence Approximations 

In this section, we will look at the approximations made in order to show classical correspondence between the 
damped PSDs of the quantum harmonic oscillator and the classical harmonic oscillator. 


1. First-Order 


We begin by using the approximations given in Eq. (81) to obtain the expression in the first line of Eq. (82) 

1 


Sxx (^) 


Tk B T 

2?71Wq 

r k B T 


2?7TWq |_ (w — wo) 2 + (r/2) 2 


1 

(w + wo) 2 + {T/2Y 

(E/2) 2 


(Jl) 


tow 2 (w 2 - w 2 ) 2 + (w - w 0 ) 2 (r/2) 2 + (w + w 0 ) 2 (r/2) 2 + (P/2) 4 


At this point, we apply the high-Q approximation for which we neglect the last terms in both the numerator and 
denominator. As well, in this approximation the only significant contributions to the PSD occur at w ss ±wo, otherwise 
the value will be small compared the the peak value. Inputting this approximation into our above PSD (we can use 
either sign, both give the same result since the classical PSD is an even function of w), we obtain 


£>xx (^0 


Tk B T 


2 Wq 


mw 2 (w 2 -w 2 ) 2 + (2w) 2 (r/2) 2 

2rfc s r 


TO [(w 2 — Wg ) 2 + w 2 r 2 ] ’ 


which is exactly the result we obtained in Eq. (68). 


(J2) 


2. Second-Order 


We follow a similar procedure here as in the previous section to show the classical correspondence of the second- 
order quantum PSD in the high-Q limit. Beginning with the high temperature approximation of the second-order 
quantum PSD given by Eq. (103), we have 


^x 2 x 2 (^) ~ 


\TOWq / 
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= rU4Y 

V. TOWg / 


(w — 2wo) 2 + r 2 (w + 2wo) 2 + r 2 uj 2 + r 2 
(w 2 + r 2 ) ^(w + 2wq) 2 + (w — 2 wq) 2 + 2P 2 ^ + ((w + 2wg) 2 + r 2 ) ((w - 2wg) 2 + T 2 ) 


(w 2 + r 2 ) ((w 2 - 4w 2 ) 2 + r 2 (w - 2w 0 ) 2 + r 2 (w + 2w 0 ) 2 + r 4 


(J3) 


As before, due to our high-Q approximation, we only concern ourselves with the PSD in the vicinity of the peaks, 
namely w « ±2wq and w « 0. As well, we neglect all terms of order T 2 or higher in the numerator and the T 4 term 
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in the denominator. With these approximations we find 

2 


64Wn 


, ) ^ p (k B T V __ 

\TOWq/ ( w 2 +r 2) (( W 2 -4 W 2) 2 + 4 r 2 w 2^ 


64r (fc B r) 2 


m 2 (w 2 + T 2 ) ((w 2 — 4wg) + 4w 2 r 2 


where again we have recovered the classical PSD that we found in Eq. (91). 


(J4) 


Appendix K: Quantum Power Spectral Density Calculations Using the Fluctuation Dissipation Theorem 


Here we calculate the first- and second-order PSDs for the quantum harmonic oscillator using the fluctuation- 
dissipation theorem, providing an independent check on our results found in Eqs. (74) and (98). 


1. First-Order 

The first-order PSD for the position operator can also be derived using the finite-temperature Green’s function 
formalism m- The main building block in this formalism is the time-ordered Green’s function for the bosonic 
operator b, 


G(r) = (T r 6(T)6t(0)>, 


(Kl) 


from which correlation functions of arbitrary operators can be constructed. Here we have introduced 7~ T , which is a 
time-ordering operator on the Matsubara contour. In particular, for the purpose of computing the first-order PSD 
for the position operator one first computes a time-ordered correlation function for x (t) m , 


n xx {r) = (TtX(t)x(O)) = £ zp f (G(r) + G(—r)), 


(K2) 


where the imaginary time r takes values between zero and the inverse temperature /? = 1/fceT. Note that we have 
also used the fact that the order of bosonic operators can be rearranged at will within a time-ordered product, and 
that for a time-independent Hamiltonian the Green’s function (Kl) depends only on the difference of its two time 


coordinates. Both the time-ordered correlation functions and the Green’s function of bosonic operators are periodic 
in imaginary time and can be given Fourier series expansions, 


G{T) = -Y J G{iu n )e~ i ^ T , 

lUJn 

H(r) = ^H(ia;„)e-^, 

iuj-n 

where ui n = 2mr/(3 , n £ Z are bosonic Matsubara frequencies [27], and the Fourier coefficients are given by 

rP 

G(iuj n ) = / G(r)e^" T dr, 

Jo 

n(iw n ) = / n(r)e laJr * T dr. 


(K3) 


(K4) 


We point out that for a single bosonic mode at frequency wo, G(iw n ) is given by 

G{iu> n ) = -—. 

CJq 

Cciii then foe determined using the fluctuution - dissipution theorem m 

S xx {u>) = 2 (n B (w) + l)Im{n^(w)}, 


(K5) 


(K6) 
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where np(w) = {e^ hu — 1) 1 is the Bose factor at frequency w and we have introduced the retarded correlation function 
II xx( u} )i defined as the analytic continuation to real frequencies of the Matsubara correlation function H xx (iiu n ), 


n ^M = dm n n (iu„), 


(K7) 


where r/ is a positive infinitesimal. To determine this retarded correlation function, we first calculate the time-ordered 
correlation function by putting Eq. (K51 into Eq. (K2) to obtain 


thus 


n xx ) ^zpf 


n ^M = a-’zpf ( 


i 


l 


iuj n + U)o iu n - W 0 


1 


1 


ui + it) + wq oj + if] — ujq 


(K8) 


(K9) 


Using the identity lim ( y + irj) 1 = V(l/y) — inS(y) for y real where V stands for the Cauchy principal value, along 

7J-V0 + 

with the relation ng(-w) = -[ng(w) + 1], the fluctuation-dissipation theorem (K6) gives 

S X x(u) = 27rx^ pf [((n) + l)<5(w - w 0 ) + (n)6(u> + w 0 )], (K10) 


where we used ?rp(wo) = (n). This result is in agreement with what we obtained for the first-order PSD in Eq. (74). 


2. Second-Order 


The second-order PSD for the position operator can also be derived in a similar way. One begins by expressing 
the position operator in terms of creation Ir and annihilation b operators. In this way, x 2 can be written as x 2 = 
xl p [(P + Q ), where we define P = bb + b' b ( and Q = 2 N + 1 with N = b' b the number operator. The operator Q 


is special in that it commutes with the Hamiltonian in Eq. (19), such that it is a conserved quantity that acquires 


no time-dependence in the Heisenberg picture, Q(t) = Q. This, in turn, implies that correlation functions of Q will 
contain a delta function at zero frequency. In terms of the operators P and Q , the real-time correlation function for 
x 2 (t) can be written as 


(x 2 (i)x 2 (0)) = xi p{ ((P(t)P( 0)) + (P(t)Q) + (QPm + (Q 2 )). (Kll) 

The expectation value of Q 2 is given by 


(Q 2 ) = ((27V + l) 2 ) = 8(n)((n) + 1) + 1, 


(K12) 


in agreement with the time-independent term in the last line of Eq. (97). The second and third terms in Eq. (Kll) 
are zero because the operator P does not conserve the phonon number N while the operator Q does. We are thus 
left with the task of evaluating the correlation function of P appearing in the expression 

(x 2 {t): r 2 (0)) = x* pi ({P(t)P{ 0)} + 8 (n)({n) + 1) + l). (K13) 

To calculate (P(t)P(O)) at finite temperature, one again defines a time-ordered correlation function, 

n pp (t) = (TrP(r)P(O)). (K14) 


The advantage of defining a time-ordered correlation function is that for Hamiltonians quadratic in the creation and 
annihilatio n op erators, as is the case here, one can use Wick’s theorem m to express npp(r) as a product of Green’s 
functions ( |Kl| ). Because the Hamiltonian commutes with the phonon number operator N, expectation values of 
products of four creation operators or four annihilation operators vanish. Using Wick’s theorem, we are left with 

Upp(t) = (r T 6(r)6(r)6t(0)6t( 0 )) + (T t P (r)P (T)b(0)b(0)) 

= 2 (G 2 (t) + G 2 (—t)), 
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where we have made the same assumptions used in determining Eq. (K21. Substituting the expressions given in 
Eq. (K3) into Eq. (K15), we obtain 


n PP {iuj n ) = - ^ ~2(G(ip n )G(iu} n - ip n ) + G(ip n )G(-iuj n - ip n )), 

Wn 


(K15) 


where p n are also bosonic Matsubara frequencies. Sums over Matsubara frequencies are most conveniently carried 
out by making use of the spectral function A{u>) 1 defined in terms of the Green’s function as 


G^u) = ^J 


A(u) 


dui. 


hence from Eq. (K5) we can read off the spectral function, 

A(u) = 2 ttS(lo — wo)- 


(K16) 


(K17) 


Substituting Eq. (K16) into Eq. (K15), we have 


n PP {iuj n ) = ^2 JI A (^)A{u}')^Yl 


1 


+ 


1 


(ip n - v)(iuJn - ipn ~ w') ( ip n - ui)(-iuj n - ip n - w 


— ) dojdoj'. (K18) 

')) 


Such sums can be performed using contour integration m, observing that the Bose factor n B (w) has single poles at 
the Matsubara frequencies u> = ip n with residue 1//3. One finds 


n PP {iu) n ) = 2 2 f (A(u])A(J)(n B {u) + n B {u') + 1) 

(27rj J J 


lUJn + UJ + UJ' lUJ n ~ UJ ~ 


- ] dui dui', 

w — ui') 


(K19) 


which upon making use of Eq. (K17| gives 


II pp{iuj n ) = 2(2 (n) + 1) 


1 


1 


iui n + 2u>o iuj n — 2wq 


Using Eq. (K71, the retarded correlation function will be given as 

tR ' ' -- • ( 1 


n $p(w) = 2(2 <n) + 1) 


+ ii y + 2wo u> + ip — 2wo / 
Substituting this result in the fluctuation-dissipation theorem ( |K6[ ) we arrive at the relation 

Spp(lu) = 47r(2(n) + 1 )(n B (ui) + 1) [<5(w — 2wo) — S(co + 2wo)]. 

Using the identity 

[»bM] 2 


we obtain 


n B (2w) = , 

2 n B {u>) + 1 

Spp(oj) = 47t [((n) + 1 ) 2 S(co — 2wq) + ( n) 2 d(ui + 2wq)] . 


(K20) 


(K21) 


(K22) 


(K23) 


(K24) 


Taking the Fourier transform of Eq. (K13) and using the definition in Eq. (|37j) , the second-order quantum PSD is 
given as 


S x 2 x 2 (u:) = x* p{ [S PP (u]) + 27t( 8(ra)((n) + 1) + l)<5(w)], 


hence using Eq. (K24| we find the second-order PSD for the position operator to be, 

S x 2 x 2 (co) = 2nx* p{ 2((n) + l) 2 <5(w — 2 wq) + 2 (n) 2 5(uj + 2wq) + (8 (n)((n) + 1) + l)(5(w) 


(K25) 


(K26) 


in agreement with Eq. (981. 
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